#### setting environment ####
require(coda)
require(ltm)

# function to compute confidence intervals
confidence <- function(x, level = 0.95, upper = TRUE){
  x <- na.omit(x)
  n <- length(x)
  mean <- mean(x, na.rm = TRUE)
  SE <- sd(x, na.rm = TRUE) / sqrt(n)
  Q <- qnorm((1 - level) / 2, lower.tail = FALSE)
  if (upper == TRUE) {
    bound <- mean + Q * SE
  } else {
    bound <- mean - Q * SE
  }
  return(bound)
}

newspaper.names <- c("Asahi", "Chugoku", "Chunichi", "Hokkaido", "Kahoku", 
                     "Mainichi", "Nikkei", "Nishinippon", "Sankei", "Yomiuri")

#### data preparation ####
# the 2017 UTokyo-Asahi voter survey data (abridged)
voter.2017 <- read.csv("voter2017_IJPP.csv")
# estimation results of the factor analysis models
load("FA_result.Rdata")

#### newspapers' ideal points ####
# posterior draws
Study2.ideal.points.draws <- t(rbind(FA.result$mcmc[[1]][, 1:10], 
                                     FA.result$mcmc[[2]][, 1:10], 
                                     FA.result$mcmc[[3]][, 1:10])) * 10
# point estimates (posterior means)
Study2.ideal.points.mean <- rowMeans(Study2.ideal.points.draws)
# 95% credible intervals (highest posterior density intervals)
Study2.ideal.points.CI <- HPDinterval(as.mcmc(t(Study2.ideal.points.draws)))

#### respondents' ideal points ####
variables <- c("Q22_1", "Q23_1", "Q23_2", "Q23_3", "Q23_4", 
               "Q23_5", "Q23_6", "Q23_7", "Q23_8", "Q23_9", 
               "Q23_10", "Q23_11", "Q23_12", "Q23_13", "Q23_14", 
               "Q23_15", "Q23_16", "Q23_17", "Q24_1", "Q24_2", 
               "Q24_3", "Q24_4", "Q24_5", "Q25_1", "Q25_2", 
               "Q25_3", "Q25_4", "Q25_5", "Q26_1", "Q27")

## handle missing data
for(i in variables) {
  n <- voter.2017[, i] == 66 | voter.2017[, i] == 99
  voter.2017[n, i] <- NA
}
voter.2017.IRT <- voter.2017[rowSums(voter.2017[, variables], na.rm = TRUE) != 0, ]

## specify variables
read.Asahi.IRT <- voter.2017.IRT$Q13_1
read.Yomiuri.IRT <- voter.2017.IRT$Q13_2
read.Mainichi.IRT <- voter.2017.IRT$Q13_3
read.Nikkei.IRT <- voter.2017.IRT$Q13_4
read.Sankei.IRT <- voter.2017.IRT$Q13_5
read.Hokkaido.IRT <- ifelse(voter.2017.IRT$Q13_1 == 99, 99, 
                            ifelse(voter.2017.IRT$PREFEC == 1 & voter.2017.IRT$Q13_8 == 1, 1, 0))
read.Kahoku.IRT <- ifelse(voter.2017.IRT$Q13_1 == 99, 99, 
                          ifelse(voter.2017.IRT$PREFEC == 4 & voter.2017.IRT$Q13_8 == 1, 1, 0))
read.Chunichi.IRT <- ifelse(voter.2017.IRT$Q13_1 == 99, 99, 
                            ifelse(voter.2017.IRT$Q13_6 == 1 | (voter.2017.IRT$PREFEC == 13 & voter.2017.IRT$Q13_8 == 1), 1, 0))
read.Chugoku.IRT <- ifelse(voter.2017.IRT$Q13_1 == 99, 99, 
                           ifelse(voter.2017.IRT$PREFEC== 34 & voter.2017.IRT$Q13_8 == 1, 1, 0))
read.Nishinippon.IRT <- voter.2017.IRT$Q13_7
not.read <- voter.2017.IRT$Q13_10

## descriptive statistics of newspaper readership (Table A.4)
national.statistics <- cbind(read.Asahi.IRT, read.Chugoku.IRT, read.Chunichi.IRT, read.Hokkaido.IRT, 
                             read.Kahoku.IRT, read.Mainichi.IRT, read.Nikkei.IRT, read.Nishinippon.IRT, 
                             read.Sankei.IRT, read.Yomiuri.IRT, not.read)
national.statistics[national.statistics == 99] <- NA

# compute the percentage of readers in prefectures where the local newspapers are distributed
local.statistics <- national.statistics
local.statistics[, c(1, 5, 6, 8, 9, 10)] <- NA
local.statistics[voter.2017.IRT$PREFEC != 34, 2] <- NA
local.statistics[! voter.2017.IRT$PREFEC %in% c(13, 16:18, 20:25), 3] <- NA
local.statistics[voter.2017.IRT$PREFEC != 1, 4] <- NA
local.statistics[voter.2017.IRT$PREFEC != 4, 5] <- NA
local.statistics[! voter.2017.IRT$PREFEC %in% 40:44, 8] <- NA

round(cbind(colSums(national.statistics, na.rm = TRUE), 
            colMeans(national.statistics, na.rm = TRUE) * 100, 
            colMeans(local.statistics, na.rm = TRUE) * 100), 1)

## estimate item parameters
IRT.result <- grm(voter.2017.IRT[, variables], IRT.param = FALSE)

# discrimination parameters (Table A.5)
discrimination.estimates <- sapply(IRT.result$coefficients, function(x) x[5])
names(discrimination.estimates) <- variables
discrimination.order <- order(abs(discrimination.estimates), decreasing = TRUE)
round(discrimination.estimates[discrimination.order], 2)

## estimate respondents' ideal points
IRT.score <- factor.scores(IRT.result, voter.2017.IRT[, variables])

## summarize the results
ideology.IRT <- data.frame(Newspaper = newspaper.names, 
                           ideology = rep(NA, 10), 
                           upper = rep(NA, 10), 
                           lower = rep(NA, 10))
# Asahi
tmp.value <- -IRT.score$score.dat$z1[read.Asahi.IRT == 1]
ideology.IRT$ideology[1] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[1] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[1] <- confidence(tmp.value, upper = FALSE)
# Chugoku
tmp.value <- -IRT.score$score.dat$z1[read.Chugoku.IRT == 1]
ideology.IRT$ideology[2] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[2] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[2] <- confidence(tmp.value, upper = FALSE)
# Chunichi
tmp.value <- -IRT.score$score.dat$z1[read.Chunichi.IRT == 1]
ideology.IRT$ideology[3] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[3] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[3] <- confidence(tmp.value, upper = FALSE)
# Hokkaido
tmp.value <- -IRT.score$score.dat$z1[read.Hokkaido.IRT == 1]
ideology.IRT$ideology[4] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[4] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[4] <- confidence(tmp.value, upper = FALSE)
# Kahoku
tmp.value <- -IRT.score$score.dat$z1[read.Kahoku.IRT == 1]
ideology.IRT$ideology[5] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[5] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[5] <- confidence(tmp.value, upper = FALSE)
# Mainichi
tmp.value <- -IRT.score$score.dat$z1[read.Mainichi.IRT == 1]
ideology.IRT$ideology[6] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[6] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[6] <- confidence(tmp.value, upper = FALSE)
# Nikkei
tmp.value <- -IRT.score$score.dat$z1[read.Nikkei.IRT == 1]
ideology.IRT$ideology[7] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[7] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[7] <- confidence(tmp.value, upper = FALSE)
# Nishinippon
tmp.value <- -IRT.score$score.dat$z1[read.Nishinippon.IRT == 1]
ideology.IRT$ideology[8] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[8] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[8] <- confidence(tmp.value, upper = FALSE)
# Sankei
tmp.value <- -IRT.score$score.dat$z1[read.Sankei.IRT == 1]
ideology.IRT$ideology[9] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[9] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[9] <- confidence(tmp.value, upper = FALSE)
# Yomiuri
tmp.value <- -IRT.score$score.dat$z1[read.Yomiuri.IRT == 1]
ideology.IRT$ideology[10] <- mean(tmp.value, na.rm = TRUE)
ideology.IRT$upper[10] <- confidence(tmp.value, upper = TRUE)
ideology.IRT$lower[10] <- confidence(tmp.value, upper = FALSE)

#### illustrating the results ####
# ascending order
ideology.IRT.order <- order(ideology.IRT$ideology, decreasing = TRUE)
# indicator for local newspapers
local <- c(0, 1, 1, 1, 1, 0, 0, 1, 0, 0)
# indicator for newspapers managing professional baseball teams
sports <- c(0, 0, 1, 0, 0, 0, 0, 0, 0, 1)

## drawing the figure
cairo_pdf("Figure_3_a.pdf", 
          width = 3, height = 3, pointsize = 7, family = "Helvetica")
par(mar = c(4, 4, 2.5, 2.5))
plot(NULL, NULL, type = "n", xlim = c(-1, 1.7), ylim = c(-0.8, 1.4), 
     main = "(a) Paper-based exposure", 
     xlab = "Newspapers' Ideal Points", ylab = "Average of Readers' Ideal Points")
segments(Study2.ideal.points.mean, ideology.IRT$lower, 
         Study2.ideal.points.mean, ideology.IRT$upper, col = "gray60", lwd = 0.5)
segments(Study2.ideal.points.CI[, 1], ideology.IRT$ideology, 
         Study2.ideal.points.CI[, 2], ideology.IRT$ideology, col = "gray60", lwd = 0.5)
points(Study2.ideal.points.mean, ideology.IRT$ideology, 
       pch = ifelse(local == 1, 23, 21), col = "gray40", 
       bg = ifelse(sports == 1, "white", "gray40"))
text(Study2.ideal.points.mean[1] + 0.11, ideology.IRT$ideology[1] + 0.05, 
     labels = newspaper.names[1], cex = 0.6)
text(Study2.ideal.points.mean[2] - 0.17, ideology.IRT$ideology[2] - 0.05, 
     labels = newspaper.names[2], cex = 0.6)
text(Study2.ideal.points.mean[3] - 0.17, ideology.IRT$ideology[3] + 0.05, 
     labels = newspaper.names[3], cex = 0.6)
text(Study2.ideal.points.mean[4] + 0.18, ideology.IRT$ideology[4] - 0.05, 
     labels = newspaper.names[4], cex = 0.6)
text(Study2.ideal.points.mean[5] + 0.15, ideology.IRT$ideology[5] + 0.05, 
     labels = newspaper.names[5], cex = 0.6)
text(Study2.ideal.points.mean[6] + 0.16, ideology.IRT$ideology[6] + 0.05, 
     labels = newspaper.names[6], cex = 0.6)
text(Study2.ideal.points.mean[7] + 0.12, ideology.IRT$ideology[7] + 0.05, 
     labels = newspaper.names[7], cex = 0.6)
text(Study2.ideal.points.mean[8] + 0.22, ideology.IRT$ideology[8] + 0.05, 
     labels = newspaper.names[8], cex = 0.6)
text(Study2.ideal.points.mean[9] + 0.13, ideology.IRT$ideology[9] + 0.05, 
     labels = newspaper.names[9], cex = 0.6)
text(Study2.ideal.points.mean[10] + 0.15, ideology.IRT$ideology[10] + 0.05, 
     labels = newspaper.names[10], cex = 0.6)
dev.off()